Real-time energy dynamics in spin- 1/2 Heisenberg chains 
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We study the real-time dynamics of the local energy density in the spin-1/2 XXZ chain starting 
from initial states with an inhomogeneous profile of bond energies. Numerical simulations of the 
dynamics of the initial states are carried out using the adaptive time-dependent density matrix 
renormalization group method. We analyze the time dependence of the spatial variance associated 
with the local energy density to classify the dynamics as either ballistic or diffusive. Our results 
are consistent with ballistic behavior both in the massless and the massive phase. We also study 
the same problem within Luttinger Liquid theory and obtain that energy wave-packets propagate 
with the sound velocity. We recover this behavior in our numerical simulations in the limit of very 
weakly perturbed initial states. 



I. INTRODUCTION 

The understanding of transport properties of low- 
dimensional systems with strong correlations still poses 
viable challenges to theorists. These include, on the one 
hand, the fundamental problem of calculating transport 
coefficients for generic models such as the Heisenberg 
chaiufii^ and on the other hand, the theoretical model- 
ing of experiments that typically require the treatment of 
spin or electronic degrees of freedom coupled to phonons, 
in particular, in the case of the thermal conductivityj^— 
Most theoretical work has focussed on the linear-response 
regime, in which the properties of current-current auto- 
correlation functions determine transport properties (see 
Refs. [l|3 for a review) . 

More recently, the out-of-equilibrium properties of one- 
dimensional systems have evolved into an active field 
of research, one reason being recent advances in experi- 
ments with ultracold atomsi^ These have paved the way 
for studying the dynamics of quantum many-body sys- 
tems that are driven far away from equilibrium in a con- 
trolled manner, with little or no coupling to external de- 
grees of freedom. Much attention has been paid to the 
question of thermalization, typically studied in so-called 
quantum quenches (see Ref. |6| and references therein). 
While global quantum quenches in homogeneous systems 



energy density <h(x,t)> 




usually do not induce any finite net currents (of either 
spin, energy, or particles), we will be particularly inter- 
ested in set-ups that feature finite net-currents. Such 
situations are realized in, for instance, the sudden ex- 
pansion of particles in optical lattices after the removal 
of trapping potcntialsi^ Further examples arc spin and/or 
particle currents induced by connecting two regions with 
opposite magnetizations or by letting two particle clouds 
collide (see, for instance, Refs. Isl-KLlI). 

Theoretical work in this context ranges from the ex- 
pansion dynamics of bosons and fermions in optical 
latticesi^"— over the dynamics of wave-packets in spin 
chainsr^r— to the demonstration of signatures of spin- 
charge separation in such set-upsj^i^ In the aforemen- 
tioned examples, non-equilibrium situations were studied 
with either finite spin or particle currents. In our work, 
we address the energy dynamics for a model that is proto- 
typical for systems with strong correlations, namely the 
spin-1/2 XXZ chain: 
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FIG. 1: (color online) Sketch of our setup: We prepare initial 
states with an inhomogeneous distribution of local energies 
and then study the time evolution of the local energy density. 



where S^ and /i = x,y,z are the components of a spin- 
1/2 operator acting on site i and S^ arc the correspond- 
ing lowering/raising operators. The global energy scale 
is set by the exchange coupling J, A is the exchange 
anisotropy in the z-direction, and L denotes the number 
of sites. Equation ([T]) describes either interacting quan- 
tum spins or, via the Jordan- Wigner transformation;^ 
spinlcss fermions. 

Specifically, wc follow the time-evolution of the local 
energy density {hi) starting from initial states that are 
far from the ground state of Eq. ([1]) and that feature 
an inhomogeneous profile in the local energy density (see 
Fig. [T] for a sketch). We emphasize that, in the main 
part of our work, we choose the initial conditions such 
that only finite energy currents exist, whereas the spin 
(particle) density is constant during the time evolution, 



hence all spin (particle) currents vanish. Obviously, an 
initial state with an inhomogeneous spin density profile 
leads to both finite spin and energy currents, and we 
revisit this case, previously studied in Refs. [20 and [2^, 
as well. 

Our work is motivated by and closely related to a spe- 
cific experiment on a spin-ladder material. Many low- 
dimensional quantum magnets are known to be very good 
thermal conductors with heat predominantly carried by 
magnetic excitations at elevated temperaturesi^i^ Ex- 
amples for materials that exhibit particularly large ther- 
mal conductivities are (Sr,La,Ca)i4Cu2404i (Refs. [3J 
and M) and SrCuOa (Ref. M)- While these experi- 
ments are carried out under steady-state conditions and 
in the regime of small external perturbations, more re- 
cently, time-resolved measurements have been performed 
on La9Ca5Cu2404i (Ref. 1371 ). For this spin ladder ma- 
terial, two approaches have been implemented: A time- 
of-flight measurement, in which one side of the sample 
is heated up with a laser pulse and the time-dependent 
response is recorded on the other side. Second, a non- 
equilibrium local heat distribution was generated in the 
surface of the material by shining laser light on it. It is 
possible to record the heat dynamics via thermal imaging 
that uses the response of an excited thin fluorescent layer 
placed on top of the spin ladder material. 

It is the latter case that we mimic in our work: 
The time evolution of local energy densities induced 
by inhomogeneous initial distributions. We utilize the 
time-dependent density matrix rcnormalization group 
(tDMRG)^"— technique. It allows us to simulate the 
dynamics of pure states whereas in the experiment, tem- 
perature likely plays a role. Our work thus addresses 
qualitative aspects in the first place, while a direct com- 
parison with experimental results is beyond the scope of 
this study: The goal is to demonstrate that in a spin- 
1/2 chain described by Eq. ([T]), the energy dynamics is 
ballistic, irrespective of how far from equilibrium the sys- 
tem is and also irrespective of the presence or absence of 
excitation gaps. To this end, we use the same approach 
as in Ref. l20|: We classify the dynamics based on the 
behavior of the spatial variance (J%{t) of the local en- 
ergy density: The ballistic case is cr|,(i) ^ t^, whereas 
diffusion implies a'^ ~ t. Our main result for the XXZ 
chain, based on numerical tDMRG simulations, is that 
energy propagates ballistically at sufficiently long times, 
independently of model parameters (such as A). One can 
then interpret the prefactor Ve in a^{t) = V^t^ as a mea- 
sure of the average velocity of excitations contributing to 
the expansion. The velocity Ve can be calculated ana- 
lytically and exactly in non-interacting models, which (in 
the absence of impurities or disorder) typically have bal- 
listic dynamics, and we consider two examples: (i) the 
noninteracting limit of the XXZ-Hamiltonian (A = 0), 
i.e., spinless fermions and (ii) the Luttinger liquid, which 
is the universal low-energy theory in the continuum limit 
of Eq. (d]) for |A| < 1. Wc show that our tDMRG rcsufis 
agree with the exactly known expansion velocity Ve in 



these two examples. 

Our main result, namely the numerical observation 
of cr|;(t) ^ t^ independently of initial conditions or 
model parameters such as the exchange anisotropy A, 
is consistent with the qualitative picture derived from 
linear-response theory. Within that theory transport 
properties of the XXZ chain have intensely been stud- 
ied in recent years, both the energy^"— and the spin 
transport 1 ^'^'^^" — Ballistic dynamics is associated with 
the existence of non-zero Drude weights. Since the to- 
tal energy current of the anisotropic spin- 1/2 chain is 
a conserved quantity for all A, the thermal conduc- 
tivity k(ll)) diverges in the zero-frequency limit and is 
given by RcK(a;) — DeS{uj), where De is the ther- 
mal Drude weight.—"— This behavior is different from 
the spin conductivity cr{ui). This quantity takes the 
form Rei7(w) = Ds5{uj) only at the noninteracting point 
A = 0, whereas for < A < 1, many numerical 
studieaii^iffii^ indicate Ds{T > 0) > 0, with a fi- 
nite weight at finite frequencies, though. Therefore, for 
< A < 1, Recr(a;) = Ds<5(a;) -I- a^c^{uj)- Recent field- 
theoretical and numerical work suggests that the reg- 
ular part (Trcg(w) of o'(cj) in massless phases is consis- 
tent with diffusive behaviori^^i^i^ A finite value of the 
current-current correlation function in the long time limit 
is associated with a finite Drude weight. Finite Drude 
weights can be traced back to the existence of conser- 
vation lawSf^i^ and in consequence, a potential rela- 
tion between integrability^ and ballistic behavior - in 
the sense of non-zero Drude weights - has been intensely 
discussed (see, e.g., Refs. Illi3 l47li60l - l62l and further ref- 
erences cited therein). Very recently, Prosen has pre- 
sented results that provide a lower bound to the spin 
Drude weight that is non-zero for A < li^ This is in 
qualitative agreement with earlier exact diagonalization 
studiesi^i^i^ The particular point A = 1 is still discussed 
controversially 42iiSi^iHi^"iSi First, no finite lower bound 
to the Drude weight is known.— and second, the qual- 
itative results of exact diagonalization studies seem to 
depend on details of the extrapolation of finite-size data 
to the thermodynamic limit and the statistical ensemble 
that is considered4ii^'^ 

Our approach that analyzes the time-dependence of 
spatial variances, albeit restricted to the analysis of den- 
sities, is numerically easily tractable and is an alterna- 
tive to the numerically cumbersome evaluation of current 
correlation functions. tDMRG has, for instance, been 
applied to evaluate current-current autocorrelation func- 
tions in the thermodynamic limit^ However, the acces- 
sible time scales are quite limited (t '^ 10/ J), making an 
unambiguous interpretation of the results difficult and 
the approach is not applicable to non equilibrium. Our 
approach allows us, at least in principle, to study the 
entire regime of weakly perturbed states to maximally 
excited ones. An earlier analysis of spin-density wave 
packets in various spin models has yielded the following 
picture (all based on the time-dependence of the spa- 
tial variance) ji22 In massless phases, ballistic dynamics is 



seen, whereas in massive ones, examples of diffusive dy- 
namics have been identified. It is important to stress that 
the observation of a variance that increases hnear in time 
is a necessary condition for the vahdity of the diffusion 
equation. 

Finahy, to complete the survey of related literature, 
recent studies have addressed steady-state spin and en- 
ergy transport in open systems coupled to baths, with no 
restriction to the linear-response regime i^MLiSiri^ These 
studies suggest spin transport to be ballistic in the gap- 
less phase of the XXZ spin chain and to be diffusive in the 
gapped phase with a negative differential conductance at 
large driving strengths. The heat current has been ad- 
dressed in Ref. IM where Fourier's law has been validated 
for the Ising model in a tilted field. 

A by-product of any tDMRG simulation is information 
on the time-evolution of the entanglement entropy. While 
this is not directly related to this article's chief case, it 
nevertheless provides valuable information on the numer- 
ical costs of tDMRG simulations. Qualitatively, speak- 
ing (see the discussion in Ref. |43 and references therein) , 
the faster the entanglement growth is, the shorter are 
the time scales that can be reached with tDMRG. Wc 
here show that the quenches studied in this work gener- 
ate a mild logarithmic increase of entanglement, which is 
why this problem is very well suited for tDMRG. Such 
a behavior is typical for so-called local qucnchesi^ This 
result might be useful for tDMRG practitioners. 

This paper is organized as follows: First, we intro- 
duce the model and the quantities used in our analysis in 
Secini Section IlII Al reviews the framework of bosoniza- 
tion, which is applied in Sec. IIIIBI to give an analytical 
derivation of ballistic spin and energy dynamics in the 
low-energy case, valid in the massless phase of Eq. ([T]). 
Sections IIVI and IVl contain our numerical results. First, 
we study the energy dynamics in the absence of spin cur- 
rents in Sec IIVI To this end, we generate an initial state 
consisting of a variable number of ferromagnetic bonds 
in the center of an antiferromagnetic chain. We calculate 
the time evolution of these states under Eq. ([I]) find- 
ing ballistic energy dynamics independent of the phase 
and the strength of the perturbation. To supplement 
these findings we derive an observable, which depends on 
the local currents, and whose expectation value is time- 
independent whenever (7^(t) '^ t^. The numerical calcu- 
lation of this quantity indicates ballistic dynamics as well. 
Section |V] revisits the scenario of Ref. I20I where local spin 
and energy currents are present during the dynamics as 
we start from states with an inhomogeneous spin density. 
In that case the energy density shows ballistic dynam- 
ics in the massless phase with a velocity matching the 
bosonization result in the limit of small perturbations. 
In the massive phase we observe a different behavior of 
the two transport channels, i.e., ballistic energy dynam- 
ics while the spin dynamics looks diffusive.— Finally, we 
summarize our findings in Sec. IVII Additionally, we dis- 
cuss the entanglement growth induced by coupling two 
regions with an opposite sign of the exchange coupling in 



the Appendix. 



II. SETUP AND DEFINITIONS 

A. Preparation of initial states and definition of 
spatial variance 

In this work, we focus on spin-1/2 XXZ chains of a 
finite length L given by Eq. ([1]) where our goal is to 
study the dynamics of an inhomogeneous distribution of 
the local energy density originating from a local quench 
of system parameters. The inhomogeneous distributions 
are generated by preparing the system in the respective 
ground states of the following Hamiltonians that are per- 
turbations of Hxxz from Eq. ([l}. First, 
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where hi is defined in Eq. ([1]), and second. 



where 



B, = Bn e ^-0 



(2) 



(3) 



(4) 



In the first case, we quench site-dependent exchange cou- 
plings. In this scenario we obtain initial states with 
large local energy densities. Typical initial states that 
are ground states of Eq. ([2]) are shown in Fig. [2j These 
states have b bonds with ferromagnetic Jt < in the cen- 
ter while the rest has antiferromagnetic Ji > 0. We refer 
to this setup as the Ji quench. 

In the second case, the dynamics is driven by an in- 
homogeneous spin density, enforced by an external mag- 
netic field applied in the initial state. This allows us to 
generate smooth spatial perturbations of (hi) with small 
differences in energy compared to the ground state of 
Eq. ([T|). We refer to this setup as the Bo quench. A 
more detailed discussion of the initial states generated 
by a Ji quench will be given in Sec. lIV Al The Bq quench 
was introduced in detail in Ref. l20l 

The definition of the local energy density from the 
Hamiltonian Eq. ([T]) is not unambiguous. For instance, 
it is always possible to add local terms to the Hamil- 
tonian whose total contribution by summation over all 
lattice sites vanishes. However, this seeming ambiguity 
can be resolved up to constants by requiring that any 
block of adjacent lattice sites X]i=; ^i i^ Hcrmitian and 
yet to have the same structure as the total Hamiltonian 
H. These details seem to be rather specific, yet for the 
definition of the appropriate local energy density within 
the Luttinger liquid description, see below, these formal 
considerations are important. For the XXZ chain the lo- 
cal energy density is therefore determined by the bond 
energies (hi). 
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FIG. 2: (color online) Profile of the local energy density (hi) in the initial states induced by a Ji quench for b = 1,3,5 [compare 
Eq. (I30|) ] for (a) A = 0.5, (b) A = 1 and (c) A = 1.5. In all cases, the system forms a region with ferromagnetic nearest-neighbor 
spin correlations in the middle of the chain. In the regions with antiferromagnetic Ji > 0, the local energy density oscillates, 
reflecting the antiferromagnetic nearest-neighbor correlations. 



To classify the dynamics of a density e^ we study its 
spatial variance. 
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<jl{t) = J2i'-fife,{t), 



(5) 



where fi is the first moment of e^ . The Bj are the normal- 
ized distribution linked to the energy density via 



e, = 5E-^{h,) 



(6) 



For diffusive behavior, we expect, from the fundamental 
solution of the diffusion equation^ that S<T^{t) ^ Dt 
grows linearly in time, where D is the diffusion constant 
(see, e.g, the discussion in Ref. [20). Within linear re- 
sponse theory the diffusion constant can be related to 
transport coefficients via Einstein relations, see e.g. Ref. 
ITOI . To be clear the observation of Sa'^ ^ t^ or 80% r^ tis 
a necessary conditions for the respective type of dynam- 
ics and time-dependent crossovers are possible. 



where {hi) = {hi) — {hi)o denotes the expectation value 
of hi in the initial state shifted by the ground state ex- 
pectation value (ft.,;)o = (^o|^i|'0o)- 
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is the energy difference between the initial state IV'init) 
[i.e., the ground state of either -ffj^jt or H^^^ ] and the 
ground state \tpo) of Eq. ^, both energies measured with 
respect to the unperturbed Hamiltonian from Eq. ([l} : 

Eq = (^'o|-ffxXz|V^o); Einit = (^init l-^XXZ l^init) • (8) 

On physical grounds, the energy density should be 
normalized by the amount of energy transported by the 
propagating perturbation. This is well approximated by 
the energy difference 6E between the initial state and the 
ground state of Eq. ([T]), as we have verified in many exam- 
ples. In some cases, though, the propagating energy is, 
on a quantitatively level, better described by estimating 
the area under the perturbations, as SE may also con- 
tain contributions from static deviations from the ground 
state bond-energies in the background. Nevertheless, 6E 
does not depend on the overall zero of energy and is an 
obvious measure of how far the system is driven away 
from the ground state. This, all together, justifies our 
definition of the e^ . 

To remove static contributions depending only on the 
initial distribution ei{t = 0), we subtract a%{t = 0) and 
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pected to grow quadratically in time in the case of ballis- 
tic behavior, where Ve has the dimensions of a velocity. 



B. Spatial variance in the non-interacting case 

For pedagogical reasons and to guide the ensuing dis- 
cussion we next calculate the spatial variance in the non- 
interacting limit of Eq. ([T|), i.e., at A = 0. Using the 
Jordan- Wigner transformation, we can write the Hamil- 
tonian as 

H = ^T.(^tS-,, + h.c.) = -lj2{4c,+, + h.c.) , (9) 

i i 

where c] creates a spinless fermion on site i. A subse- 
quent Fourier transformation diagonalizes the Hamilto- 
nian: 
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Since we will compare with numerical results on systems 
with open boundary conditions, we obtain 



Efc ~ — Jcos(fc); k 
Next, we compute 
5al{t)=Y.^,{t){^^ 
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with Bi from Eq. © and hi ~ — J(c,]ci_|-i + h.c.)/2. By 
expressing c\ through their Fourier transform and by 



plugging in the time evolution of c^, 
after straightforward calculations: 
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we finally obtain. 
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i.e., ballistic dynamics independently of the initial state. 
Terms linear in t will be absent if in the initial state, the 
density is symmetric with respect to its first moment, i.e., 
s^+5 = e^_5 and if the wave packet has no finite center- 
of-mass momentum at t = already. In the remainder 
of the paper we will work under these two additional 
assumptions that are valid for all initial states considered 
in our work. The prefactor Vj is given by: 



^E = J^YI '^kvlSuk 



(13) 



C. Energy current 

Another aspect worth noting is that the time-evolving 
state carries a nonzero energy current, a situation that 
usually does not appear in the case of a global quench. 
From the equation of continuity for the energy density 
one can derive the well-known expression for the local 
energy current operator— 



jf = J2 5,_i ■ (5, X 5,+i) , 



(18) 



where Vk — dek/dk and 
5nk 



^init 



nk 



is the difference between the momentum distribution 
function (MDF) in the initial state and the one in the 
ground state of Eq ([1]). Since we use open boundary 
conditions, we compute Uk from 
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— \^sin (/cr) sin (fcr')(cJ.Cr') . (14) 



We can also express 5E via Suk'- 

5E = '^ek Suk 



The expression Eq. ([T3| suggests that Ve is the average 
velocity of excitations contributing to the propagation of 
the wave packet. Characteristic for ballistic dynamics, 
yj is fully determined by the initial conditions through 
5nk- 

For completeness, we mention that an analogous cal- 
culation can be done for the spatial variance as of the 
spin density. This quantity is defined as 



^^W^=77E(*-/^)'-(^^'W + l/2). 



(15) 



The normalization constant TV measures the number of 
propagating particles. The spin density is, in terms of 
spinless fermions: 
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c\c^ - 1/2 = n,- 1/2 . 



The result for the spatial variance of the spin density is 
5al{t) = al{t)-al{Q) = Vlt' (16) 



with 
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(17) 



Although we started from the Hamiltonian for A = 0, 
we stress that Eqs. ([12]), dUl), (HS}, and (HH) are vahd for 
any dispersion relation e^., irrespective of the presence of 
a gap, provided that k has the meaning of a momentum. 



where S = {S^, S'^ ,AS^). With periodic boundary con- 
ditions, the total current Je = T^ .- 7,^ is a conserved 
quantity, i.e., [H,Je] = (see Ref.l43l). On a system 
with open boundary conditions such as the ones that 
are well-suited for DMRG, this property is lost, yet the 
dynamical conductivity still has a quasi-Drudc peak at 
very low frequencies, reminiscent of the true Drudc peak 
ReK{u) = De5{uj) of a system with periodic boundary 
conditions.— The latter form is recovered on a system 
with open boundary conditions as L ^ oo (Ref. 1681 ). 
showing that ballistic dynamics due to the existence of 
globally conserved currents can still be probed on sys- 
tems with open boundary conditions. 

To connect the local energy currents to the spatial vari- 
ance of the time-dependent density one can rewrite the 
time derivative oi a'^it) using the equation of continuity, 
assuming no current flow to sites at the boundary (this 
assumption is justified in our examples as long as we re- 
strict ourselves to times before reflections occur at the 
boundary in our simulations): 
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= -0f)+E(2'--2A^ + l)0V''(i))- (19) 

If cr|(t) = V|t2 + h and fi ^ /i(t), then using (Je) = 
leads to: 

^ 1 

Y, r dt{j^{t)) ^ -d?a%{t) = y| = const . (20) 



If we interpret this equation as an operator equation, 
then we see that we can define a quantity J^ via: 
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If for a given initial state and over a certain time window, 
(J|;(t)) = const, then we have identified a regime with 
ballistic dynamics, S(j'^{t) ^ P. If {JE{t)) =const holds 
for all times and initial states, then J^ is a conserved 
quantity, [H, J'^] = 0. This is the case at A = 0, the 
non-interacting limit of Eq. ([ij, where {J'^) = V^SE 
from Eq. P^ . 



We emphasize that we have here identified a opera- 
tor that connects the phenomenological observation of a 
quadratic increase of a'^it) to the local energy currents. 
In ballistic regimes, its expectation value becomes sta- 
tionary. 

For completeness, we mention an analogous result in 
the diffusive regime where a'^ ^ t. Then, expectation 
values of the operator 



JE = J2('-^^)Jrit) 



(22) 



are time independent. Obviously, similar expressions can 
be written down for the spatial variance associated with 
the spin density. 



III. PROPAGATING ENERGY AND SPIN 
WAVE-PACKETS IN A LUTTINGER LIQUID 

In the gapless phase, i.e., for |A| < 1, the low-energy 
and low-momentum properties of the XXZ chain can be 
described by an effective Luttinger liquid theory;^ In the 
following we want to analyze the energy density and the 
spin dynamics of the XXZ chain in this exactly solvable 
hydrodynamic limit. Specifically, we show that at least 
asymptotically for large times, the spatial variance al- 
ways grows quadratically both in the case of spin and 
energy dynamics. In addition, we work out the precise 
dependence of the prcfactor in front of the t^ increase 
of the spatial variance on system parameters. Since our 
DMRG results to be presented in Sees. |IV] and |V] show 
that (j'^{t) ~ t^ at any A, we did not investigate the in- 
fluence of marginally relevant perturbations at A = 1 on 
the wave-packet dynamics. In passing, we mention that 
in the massive phase, where the appropriate low-energy 
theory is the sine-Gordon model, the expansion velocity 
could also be derived at the Luther-Emery point (this 
case was studied, in, e.g. Refs. [25||2£ 



A. Bosonization of the anisotropic spin-1/2 chain 

The Hamiltonian Eq. ([l]) can be mapped onto a system 
of interacting spinlcss fermions via the Jordan- Wigner 
transformation!^ Within a hydrodynamic description in 
terms of a linearized fcrmionic dispersion relation the 
Hamiltonian can be represented in terms of a Luttinger 
liquid theory (LL) 
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(23) 



using the notation of Ref. \T^ The sum of the two left- 
and right-mover densities pl{x) -\- pb,{^) of the spinless 
Jordan- Wigner fermions is proportional to the contin- 
uum approximation of the local magnetization Sf up to 
a constant. The sound velocity u can be related to the 



parameters of the XXZ chain in Eq. ([T]) via the group 
velocity2^ 



J 



TT sin(i/) 



(24) 



with cos J/ — A. Similarly, the Luttinger parameter K is 
given by the relation K = 'k/\2{\ ~ v)]. In the noninter- 
acting case A = we have K ^1 and u = J . 



B. Ballistic dynamics in the gapless phase 

Within the Luttinger liquid description for A < 1, 
an initially inhomogeneous local energy density profile 
always propagates ballistically independently of the de- 
tails of the perturbation as can be seen from general ar- 
guments. For the effective low-energy Hamiltonian the 
probability distribution e(x,t) associated with the local 
energy density is given by 



e{x,t)=£ ^{'>Jjinit\h[x,t)\'4!init) 
where |'0init) is the initial state, 

h{x) = u{K + X-i)/(87r) ^ d^,^l{x)d^^^{x) 



(25) 



and 
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Jt) 



(26) 



(27) 



For the exact definition of the fields lys},, see, e.g. 
Ref. [72|. The local energy density operator consists of 
decoupled left- and right-moving contributions in the ba- 
sis in which the Hamiltonian for the time evolution is 
diagonal. This allows for a separation of e(x, t) into 
left- and right-moving contributions which both propa- 
gate with the sound velocity Vg-. e{x, t) = eL{x + Vgt, t = 
0) + eR{x-Vgt,t = 0). 

Assuming a L -^ R symmetry in the initial state, i.e., 
a state with zero total momentum, one obtains for the 
variance from Eq. ([5]) 



6al{t) = a%it)-aUt = 0) = {VEt)' 



(28) 



for all times t with Ve = Vg ^ u. This results can also 
be obtained from evaluating Eq. p3p in the continuum 
limit. 

In the case of an initial L <r^ R asymmetry in the ini- 
tial state we get Sa^{t) — > (vgt)'^ for t ^- oo, but the 
short time behavior may differ. Thus, within the validity 
of a Luttinger liquid description the energy transport is 
always ballistic for all initial conditions. This is evident 
from a physical point of view as all excitations propagate 
with exactly the same velocity Vg, the left- movers to the 




FIG. 3: (color online) Energy difference SE between the initial 
state and the ground state for the Ji quench as a function of 
b for A = 0.5, 1, 1.5. The inset shows the hierarchy of states 
with increasing total spin which appear as initial states when 
the total spin is a good quantum number, i.e., at A = 1. 



left and the right-movers to the right. Note that the ap- 
plicability of a Luttinger liquid description is manifestly 
restricted to cases in which the initial energy density pro- 
file is a smooth one in the sense that the associated exci- 
tations do not feel the nonlinearity of the fermionic dis- 
persion relation. Thus, the time-evolution starting from 
initial profiles such as the ones shown in Fig.[5]are beyond 
the scope of this low-energy theory. 

In analogy to the above arguments, the dynamics of 
spin-density wave-packets is also ballistic in the XX Z 
chain for A < 1 in the Luttinger liquid limit. In the 
bosonic theory, the spin density is proportional to pl{x)+ 
Pr{x) up to a constant, see Sec. IIIIAI The associated 
probability distribution p{x,t) = Q~^{pl + Pr)/'^t^i with 
Q ~ J dx{pL + pr) /2Tr, can again be separated into a 
left- and a right- moving contribution, i.e., 

p{x, t) ^pl(x + Vgt, t = 0) + pr{x - Vgt, < = 0) . (29) 

Thus, similar to the case of the energy dynamics, one 
finds ballistic behavior for |A| < 1 consistent with the 
numerical results of Ref. uO- 



IV. DMRG RESULTS FOR THE Jt QUENCH 

Now we turn to the numerical simulations. Using the 
adaptive time-dependent DMRG^"— method we can ac- 
cess the real-time dynamics of initial bond energy distri- 
butions. Within this approach we can probe the micro- 
scopic dynamics including the time dependence of bond 
energies or the entanglement entropy starting from vari- 
ous initial states in an essentially exact manner without 
limitations in the range of parameters. We discuss the 
pure energy dynamics in the absence of spin currents in- 
duced by the Jt quench in this section. We detail the 
construction of initial states and their specific features, 
then move on to the analysis of the time evolution of the 
bond energies. We calculate the spatial variance and the 
related quantity J'^ and discuss the emergent velocities 
of the energy dynamics. Within the numerical accuracy 



of our simulations we find a quadratic increase of a^^ (t) in 
all cases studied. However, it seems that for a Ji quench a 
large number of different velocities contribute as opposed 
to the Luttinger Liquid theory result, the latter valid at 
low energies. Our study of the energy current during the 
time evolution and the time evolution of the expectation 
value (J|;(i)), defined in Eq. (|2T|) . gives additional in- 
sights into short-time dynamics and further validates the 
conclusion of ballistic energy dynamics. 



A. Initial states 

Let us first describe the typical shape of initial states 
induced by a Ji quench on a few bonds in the middle of 
the spin chain. To be specific, in the Hamiltonian Eq. ([2]) 
we set 



J i<L/2-b 

J,= { -J for L/2-b<i< L/2 + b 

J i>L/2 + b 



(30) 



which provides us with initial states with an inhomoge- 
neous energy density profile with a width of 2b of the fer- 
romagnetic region. Outside this ferromagnetic region we 
obtain antiferromagnetic nearest-neighbor correlations. 

Figure [2] shows the profile of the local energy density 
of XXZ-chains with L = 100 sites with (a) A = 0.5, 
(b) A = 1 and (c) A = 1.5, induced by a sign change 
of Ji on b = 1,3,5 bonds [compare Eq. ([30])]. obtained 
using DMRG with m = 200 states exploiting the C/(l) 
symmetry to ensure zero global magnetization S^^^ = 
^j(5'f) = and, in consequence, (S'f) = 0. In all cases 
shown in Fig. [21 the system forms a region with ferro- 
magnetic nearest-neighbor correlations in the middle of 
the chain. Note that for A 7^ 1, (hi) is the sum of the 
nearest-neighbor transverse and longitudinal spin corre- 
lations, the latter weighted with A. In the regions with 
antiferromagnetic Ji > 0, the local energy density oscil- 
lates, refiecting the antiferromagnetic nearest-neighbor 
correlations. Figure [3] shows the energy difference SE. 
As a function of 6, the energy difference SE increases 
linearly once the smallest possible ferromagnetic region 
has been established. The minimum energy difference 
SE — ii^init ~ Eq is of thc order of 2 J, i.e., initial states 
that are only weak perturbations of the respective ground 
state cannot be generated using a Ji quench. 

At the isotropic point A = 1, we can explain the de- 
pendence of the initial state on the width 6 in a trans- 
parent manner. Thc ground state energy per site for the 
antiferromagnetic ground state is known from the Bethe 
Ansatz to be hm^^^oo Eo{L)/L = -hi(2) + 1/4^ while 
for the ferromagnetic ground state Eq/{L — 1) = 1/4 
excluding the boundary sites which gives rise to a very 
small system-size dependence. By growing the ferromag- 
netic region symmetrically with respect to the center of 
the chain and taking Eq{L) from the unperturbed ground 
state with open boundaries, we obtain states with an en- 
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FIG. 4: Time evolution of the bond energy distribution starting from initial states with 6=1 from Fig. [2] for (a) A — 0.5, (b) 
A = 1 and (c) A = 1.5. Despite the different ground state phases, for the selected values of the exchange anisotropy A, main 
features of the dynamics such as two distinct rays extending from the edges of the perturbation are similar. The solid white 
lines for A = 0.5 and A = 1 indicate the propagation of a single excitations starting in the middle of the chain at time i = 
moving with the group velocity Vg from Eq. (|24|) . This is also the velocity in the outer rays. 



ergy that increases as 

SE{b) = (2*6-1)- (Eo/iL - 1) - 0.25) + SE° , (31) 

for our finite system size {SE*^ is simply an ofF-sct). Equa- 
tion ([31]) exactly reproduces the data for A = 1 shown 
in Fig. 13 

Furthermore, at A = 1, the total spin 



o2 



E^^-E^. 



(32) 



is a conserved quantity. Since the ground state calcu- 
lation only respects the conservation of magnetization 
(•^tot = 0) ■^•^ obtain a hierarchy of states with S > 0. 
This can be easily understood by considering the block 
structure of the initial state. Taking, e.g., a total of 
i = 100 spins and assuming a ferromagnetic region of 
only two spins (i.e., 6=1), the two ferromagnetic spins 
are fully polarized with a total spin of 5 = 1 while each 
of the antiferromagnetic blocks have 49 spins and there- 
fore a total spin of S* = 1/2. Thus, the total spin of the 
whole chain is 5'tot = 2. Increasing the width of the fer- 
romagnetic region by one, i.e., to b = 2, we have S = 2 in 
the middle and the antiferromagnetic blocks are of even 
length, which both have 5* = in their ground state. 
This pattern repeats itself upon increasing the length 26 
of the ferromagnetic region. 



B. Time evolution of bond-energies after a J; 
quench 

Now we focus on the time-evolution of the local energy 
density induced by the aforementioned perturbation. At 
time i = 0+, we set all J,; = J and then evolve under the 
dynamics of Eq. ([1]). The DMRG simulations are car- 
ried out using a Krylov-space based algorith m^^i^^ with 
a time-step of typically 0.25J and by enforcing a fixed 
discarded weight. We restrict the discussion to times 
smaller than the time needed for the fastest excitation to 
reach the boundary. 



1. Ji quench: Qualitative features 

Figure [4] shows the time evolution of the bond energies 
{hi{t)) as a contour plot for A = 0.5, 1, 1.5 at 6 = 1. De- 
spite the different ground states for the selected values 
of anisotropy, all features of the dynamics such as two 
distinct rays starting at the edges of the block of ferro- 
magnetic correlations, are similar. The solid white lines 
for A = 0.5 and A = 1 indicate an excitation spreading 
out from the center of the ferromagnetic region with the 
group velocity given by Eq. ([M|) (these lines are paral- 
lel to the outer rays visible in the figure, i.e., the fastest 
propagating particles). Note that Eq. ([M]) holds only in 
the gapless phase (|A| < 1). Besides the outer rays that 
define a light cone structrue. Fig. 2] unveils the presence 
of more such rays inside the light-cone. Since our partic- 
ular initial states have a sharp edge in real space, there 
ought to be many exciations with different momenta k 
contributing to the expansion. 



2. Ji quench: Spatial variance 

Our main evidence for ballistic dynamics in both 
phases is based on the analysis of the spatial variance, 
shown in Fig. [5j Fitting a power-law (straight lines) to 
the data, i.e., (J^{t) — cr^{0) = at^ yields a quadratic 
increase with /3 ~ 2, classifying the dynamics as ballistic. 

In order to estimate uncertainties in the fitting param- 
eter a, we compare this to the results of fitting a pure 
parabola a% [t) — a% (0) = V| i^ to the data. Typically, 
yj deviates from a by about 10% while the exponent 
of the power-law fit is usually different from 2 by 5%. 
As an example, for A = 0.5, 6 = 1 we obtain /3 = 2.03 
and a = 0.53 vs. Vj = 0.6J^. The main reason for the 
deviation of /3 from 2 is, in fact, that the short time dy- 
namics is not well described by a power law at all over a 6- 
dependent time-window. We shall see later, in Sec. IIV C[ 
that the ballistic dynamics sets in only after the block 
of fcrromagnctically correlated bonds has fully 'melted'. 
Indeed, by excluding several time steps at the beginning 




FIG. 5: (color online) Spatial variance of the evolving energy 
distribution for (a) b = 1, (b) 6 = 5 and A = 0.5, 1, 1.5. 
Fitting a power- law (straight lines) to a%{t) — cte{0) = ctr 
yields a quadratic increase with sufficient accuracy, classifying 
the dynamics as ballistic. For instance we find a = 0.53, /3 = 
2.03 for A — 0.5 and 6=1 [black circles in (a)]. We do not 
find any qualitative difference between the massless (|A| < 1) 
and the massive (A > 1) phase. The deviations between the 
fit and the tDMRG data in the A = 1.5 curves at the largest 
times simulated are due to boundaries. 
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FIG. 6: (color online) Long-time evolution exploiting the con- 
servation of total spin Stot at A = 1 for L = 200 sites using 
an initial state with 6=1. For comparison we plot the result 
for L = 100 sites using only f/(l) symmetry. Fixing the dis- 
carded weight to 10~^ we need less than half the number of 
states. Furthermore, we find that the spatial variance is very 
robust against finite-size effects. 



time evolution respecting SU{2) symmetry (blue trian- 
gles) for a system of L = 200 sites and A = 1, 6 = 1 
compared to the result from Fig. [5] for L = 100 sites 
(red squares). We still find a quadratic increase of fT|,(t) 
and thus ballistic dynamics for times up to t ~ 60/J 
and in addition, the prcfactor does not depend on the 
system size. Both simulations were carried out keeping 
the discarded weight below 10~^ which requires at most 
m = 900 states using only C/(l) symmetry on L = 100 
sites versus a maximum oi m = 400 using SU(2) for 
L = 200 sites. 



of the evolution from the power-law fit, we observe that 
/3 ^- 2 and a ^- Vj. Therefore, we will present results 
for V^, obtained by fitting a%{t) - cr%{0) — V^t"^ to our 
tDMRG data. 



3. Exploiting SU{2) symmetry ai A = 1 for the Ji quench 

Before proceeding to the discussion of the expansion 
velocity Vj, we wish to discuss the long-time limit, which 
can be accessed in the case of A = 1. Since our per- 
turbation is proportional to the operators for the local 
energy density, global symmetries of the unperturbed 
Hamiltonian are respected by the initial states of the type 
Eq. (PO]) . Therefore, at A = 1, we can exploit the con- 
servation of total spin S, a non-Abelian symmetry. This 
can be used to push the simulations to much longer times, 
since we can perform the time evolution in an SU{2) in- 
variant basis.— The number of states needed to ensure 
a given accuracy is reduced substantially compared to a 
simulation that only respects U{1) symmetry. Therefore, 
we can work with larger system sizes and study the long- 
time dynamics of the energy density. As we can reach 
longer times, we can also analyze and discuss finite-size 
effects for A = 1 here. Figure H] shows our result for the 



4-. Expansion velocity 

The results for V^ are collected in Fig. [7] and plotted 
as a function of 6E for A = 0,0.5,1,1.5. In the non- 
interacting case, A = 0, y| is constant for b > 2, while 
at & = 1 (the smallest possible SE), V^ = 0.5J^. For all 
A > 0, yj slightly decreases with SE and V| is much 
smaller than v^ given by Eq. ([M|) . suggesting that indeed, 
many velocities contribute during the expansion of the 
energy wave-packet. 

Intuitively, one might associate the decrease of I/J, 
which is a measure of the average velocity of propagating 
excitations contributing to the expansion, to band curva- 
ture: The higher SE, the more excitations with velocities 
smaller than Vg are expected to factor in. 

It is instructive to consider the non-interacting limit 
first by comparing the numerical results obtained from a 
time-evolution with exact diagonalization to the analyt- 
ical (and also exact result) from Eq. ([T^ . To that end 
we need to compute the MDF [see Eq. p^ ] of the initial 
state. Our results for A = 0, which are shown in Fig. |51 
unveil a peculiar property: The Ji quench always induces 
changes at all k, i.e., the system is not just weakly per- 
turbed in the vicinity of kp. This is not surprising since 
our initial states have sharp edges in real-space (compare 
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FIG. 7: (color online) Prefactors V| of the fits ct| - cr|(0) = 
Vji'^ as a function of SE for A = 0, 0.5, 1, 1.5 and Ji quenches 
with b = 1,2,3,4,5 (for A = 1.5, we show b = 1,2,3 only). 
For A > 0, Vl decreases slightly with b while V| < Ug. At 
A = 0, Vj is roughly constant for 6 > 2. 



FIG. 8: (color online) MDF of the initial states generated by 
a Ji quench at A = 0.5 and A = (inset), with b — 1,3,5. 
For comparison we include the MDF of the groundstate (solid 
black line). 



Fig. [5]). Moreover, the Ji quench changes the MDF in 



such a way that Snk{b) 



-.init 



(b) — rifc is point-symmetric 



with respect to kp = tt/2 , where kp is the Fermi wave- 
vector. As Fig.[7]shows, V| as extracted from fits to Sa^ 
(solid symbols) and Vp from Eq. (|13p [open symbols] per- 
fectly agree with each other, as expected. 

The MDF of initial states for the interacting systems 
are also such that Suk ^ at all momenta and we may 
therefore conclude that the observation Vp < Vg is due 
to the fact that the Ji quench induces many excitations 
with velocities smaller than Vg (compare the data shown 
for A = 0.5 shown in Fig. [8]). Of course, Eq. (fT3)) is not 
directly applicable to the interacting case since, first, it 
does not account for the correct eigenstates at A 7^ and 
second, in general, (hi) ^ {J{S^S~^^i +h.c.)/2). Never- 
theless, by numerically calculating Suk for the interacting 
system and by using the renormalized velocity in Eq. (J13p 
instead of J [i.e., J -^ Vg{A)], we obtain an estimate for 
Vp from 



V^ 



a 
SE 



y^ cos(fc) sin^ {k)Snk 



(33) 



This reproduces the qualitative trend of the tDMRG re- 
sults for Vp as we exemplify for A = 0.5 in Fig. [71 

To summarize, the overall picture for the time evolu- 
tion of the bond energies after a Ji quench is: Energy 
propagates ballistically with an expansion velocity Vp 
that is approximately given by Eq. p3p . Combined with 
the observation that on a finite system, a Ji quench in- 
duces changes in the MDF at all momenta k, we conclude 
that many excitations contribute to the wave-packet dy- 
namics, resulting in Vp < Vg, both in the non- interacting 
and in the interacting case. 



C. Energy currents 

To conclude the discussion of the Ji quenches we 
present our results for the local energy currents at A = 1 



in Fig. ini By comparison with Fig. Hl^b), we see that 
the local current is the strongest in the vicinity of the 
wave packet. The energy current in each half of the 
system becomes a constant after a few time steps, i.e., 
j|'/2 := J2i=i~ jf reaches a constant value. We plot 
the absolute value of (^f/2) for A = 0.5, 1, 1.5 for 6 = 1 
in Fig. [TUT a). The qualitative behavior is independent of 
A: As soon as the initial perturbation has split up into 
two wave packets, we have prepared each half of the chain 
in a state with a constant, global current (j£',2)=const. 
For a system with periodic boundary conditions, the total 
current Je = 'Ylii iF i^ ^ conserved quantity4^ Since the 
effect of boundaries only factors in once these are reached 
by the fastest excitations, we directly probe the conser- 
vation of a global current with our set-up, after some ini- 
tial transient dynamics. Therefore, we can link the phe- 
nomenological observation of ballistic wave-packet dy- 
namics to the existence of a conservation law in the sys- 
tem. 

While the currents {Jp/2) clearly undergo some tran- 
sient dynamics [see Fig. [TUTa)]. we have derived a quan- 
tity in Secini called J^,, whose expectation value is sta- 
tionary if cr|, ^ t^ . Wc now numerically evaluate {J pit)) 
from Eq. (|2ip which provides an independent probe of 
ballistic dynamics. Figure rTOT b) shows our results for 
A = 1 and Ji quenches with 5=1,2, 3, 4, 5. It turns out 
that {Jp{t)) is indeed constant at sufficiently large times. 



t^. InSecHVH 



consistent with the observation of 6a p 
we have noted that Sa^ / t^ at short times t < b/J. 
This renders {Jp(t)) a time-dependent quantity over the 
same time window: Clearly, the time window over which 
{Jp{t)) 7^ const, depends on b [see Fig.[TU|b)], which sug- 
gests that the deviation of ballistic dynamics is associated 
to the 'melting' process of the region with ferromagnetic 
correlations. We have carefully checked that these ob- 
servations are robust against errors in the calculation of 
time derivatives in Eq. ([2T|) induced by the finite time 
step. Since {Jp{t)) is time-dependent (at least at short 
times), we conclude that J^ is not a conserved quantity 
in the interacting case. Finally, within our numerical ac- 
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FIG. 9: (color online) Real-time evolution of the local energy 
current Eq. p8|l at A = 1 for a Ji quench with b — 1. 



curacy and as an additional consistency check, we find 
that {J'^)/5E = a in the stationary state as expected 
from the discussion in Sec. Ill CI 

To summarize, (J|;(i))=const whenever 5a'^ ~ t^ but 
( J|;) is very sensitive to the initial transient dynamics in 
the energy dynamics and becomes constant after a time 
K, bJ. Furthermore, our set-up serves to prepare each 
half of the system in a state with a finite global energy 
current {Jfi^ that, after some transient dynamics, does 
not decay since the global energy current operator is a 
conserved quantity. 



COUPLED SPIN AND ENERGY DYNAMICS 




FIG. 10: (color online) (a) Absolute value of the current in 
each half of the system. A constant value is reached after 
t f» 5/ J. (b) The quantity {J%{t)) from Eq. ([21]) derived from 
a pure quadratic increase of the spatial variance for A = 1 and 
b = 1,2,3,4,5. This quantity is constant, as expected from 
the discussion in Sec. IIVCI except for the initial transient 
dynamics at t < b/J. 



After focussing on the energy dynamics in the absence 
of spin- /particle currents we now revisit the case of spin 
dynamics starting from states with {S^{t = 0)) y^ 0. 
Thus, during the time evolution, the local spin and en- 
ergy currents arc both non-zero. In Rcf. l20l . the dynam- 
ics of the magnetization was studied, where the inhomo- 
geneous spin density profile was induced by a Gaussian 
magnetic field in the initial state. We take the initial 
state to be the ground state of Eq. ^ in the sector with 
zero global magnetization, i.e., S^^^ = J^ii'^i) = 0- Such 
a perturbation naturally also results in an inhomogeneous 
energy density in the initial state, which is coupled to the 
spin dynamics during the time evolutioui^ 



A. Massless phase 

In Fig. Illf a) we compare the initial magnetization 
(black solid line) and the local bond energies (dashed red 
line) induced by a Gaussian magnetic field with Bq ^ J 
and (To = 5 a-t A = 0.5 finding qualitatively the same 
pattern: Both the spin and the energy density follow the 
shape of the magnetic field, resulting in a smooth per- 
turbation with small oscillations in the background away 
from the wave packet. 

For the time evolution of the bond energies at < A < 
1, we perform an analysis of their spatial variance analo- 



gous to the discussion of the Ji quench, finding ballistic 
dynamics in the massless phase. Since with a Bq quench, 
initial states with very small SE can be produced, we 
next connect our numerical results to the predictions of 
LL theory, valid in the limit 5E <^ J (compare Sec. lIIII) . 

Since we enforce zero global magnetization, we draw 
magnetization from the background into the peak.— 
Therefore, one has to carefully estimate the contributions 
to SE that do not contribute to the time-dependence of 
bond energies yet change the background density Ubg. 
The latter, in turn, affects the expected group veloc- 
ity and we thus expect to recover the LL result derived 
for the half-filled case, i.e., propagation with Vg from 
Eq. p4)). in the limit of large systems where n^g — )• 
1/2. Furthermore, Bq quenches induce 2fci^-oscillations 
in the spin and energy-density.— To account for this we 
use coarse graining, i.e., averaging the energy density 
over neighboring sites, and we take the sum only over 
the area of the peak when estimating SE: We obtain 
^^pcak _ J2l^^l+l{{h,) - {h,)o) where {h,)o denotes the 
ground state expectation value. From this quantity we 
calculate the velocity via V| — > V| ■ SE/SE^°^^, which 
is shown in Fig. [TlT b). Note that while SE-p^ak is the 
correct normalization to obtain the correct velocities, we 
label our initial states via SE. At A = (blue circles), Fj 
decreases linearly as a function of SE. Next we compare 
the result from the low-energy theory from Sec. lIIII (solid 
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FIG. 11: (color online) (a) Magnetization (solid black line) 
and energy density (dashed red line) in the initial state, for a 
-Bo quench with Bq = J and ao = 5 for A = 0.5 on a lattice 
of L = 200 sites, (b) Prefactor Vj of <5cr|(t) = V|t^ for the 
energy dynamics after a Bo quench in the massless phase of 
the XXZ chain, compared to the group velocity [Eq. (|24|) ] for 
A = 0, 0.5 and L — 200. On this system size and in the 
limit of small perturbations, Vj is approximately 5% smaller 
than the prediction from the Luttinger Liquid theory for both 
A. For A = 0, finite-size scaling of Vi(SE —>■ 0) using L — 
100, 200, ..., 800 yields V| -^ vl as shown in the inset. 



symbols at 5E = 0) to our tDMRG data. For both A = 
and A = 0.5, Vj for L — 200 sites is approximately 5% 
smaller than v'i from Eq. ([M]) . which is mainly due to 
the deviation of the background density from half fill- 
ing. While it is hard to get results for larger systems 
than L ~ 200 in the interacting case, we can solve the 
A = case numerically exactly in terms of free spinless 
fermions, allowing us to go to sufficiently large L to ob- 
serve V|(L) -^ v^ as L ^ oo. The inset of Fig. [TlT b) 



shows the finite-size scaling of V^{L) for A = using 
L = 100,200, ...,800 which yields V| -^ v^ in the limit 
L — > cx), taking first SE -^ for each system size. We 
thus, in principle, have numerical access to the dynam- 
ics in the low energy limit well described by Luttinger 
Liquid theory using a Bq quench. 



B. Massive phase 

In Ref. [23 examples of a linear increase of the spatial 
variance of the magnetization cr|(i), defined in Eq. ((TBI) , 
were found in the massive phase, which were interpreted 
as an indication of diffusive dynamics. We now demon- 
strate that while the spin dynamics may behave diffu- 
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FIG. 12: (color online) Time dependent bond energies for the 
dynamics induced by a Bo quench with Bo = 1.5J, ao = 5 
on a chain of L = 200 sites at A = 1.5: In this case, both 
local spin and local energy densities are perturbed and the 
corresponding local currents are non-zero. 



sively, i.e., Sag ^ t over a certain time window, the en- 
ergy dynamics in the same quench is still ballistic, i.e., 



5^1 



t^ 



In Fig. [T21 we show the full time evolution of the bond 
energies for a Gaussian magnetic field with Bq — 1.5 J 
and ctq = 5 on a chain of L = 200 sites at A = 1.5. It 
consists of two rays propagating with opposite velocities. 
In Fig. [131 we compare the spatial variance of the mag- 
netization cr|(i) to the one of the bond energies cr'^ii) 
calculated in the same time evolved state. The main 
panel of Fig. [T51 shows <J^^{t) — cr|;(0) which is very well 
described by a power- law fit with an exponent /3 = 1.98 
on the accessible time scales. The inset of Fig.[T51displays 
the data for ^cr|(i) = cr|(i) — cr|(0) taken from Ref. [20. 
The spatial variance of the energy density is quadratic 
in time, even at times t > 12/J where the spatial vari- 
ance of the magnetization increases only linearly. This 
example reflects the qualitative difference between spin 
and energy transport in the massive phase of the XXZ 
model at zero global magnetization: The conservation of 
the global energy current is consistent with the obser- 
vation of ballistically propagating energy wave-packets 
while spin clearly does not propagate ballistically. Our 
result, obtained in the non-equilibrium case with a zero- 
temperature background density, is consistent with the 
picture established from both linear-response theor y^^i^^ 
and steady-state simulations <22^i^ 

Very recently, Jesenko and Znidaric have also studied 
the time evolution of spin and energy densities induced 
by a Bq quenchi^ They concentrate their analysis on 
the velocity of the fastest wave-fronts, contrasting energy 
against spin dynamics. Based on the presence of these 
rays of fast propagating particles, they claim that the 
wave-packet dynamics still has ballistic features. How- 
ever, their analysis neglects the influence of slower exci- 
tations that also contribute to the dynamics of the wave 
packet, which is captured by the variance, and it ignores 
the decay of the intensity in the outer rays that we typi- 
cally observe whenever Sog ~ ti^ The latter is, if at all, 
weak in a ballistic expansion characterized by (5cr| ~ t^ . 
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FIG. 13: (color online) Spatial variance of the energy density 
(main panel) and the spin density (inset), induced by a Bo 
quench with Bo/ J — 1.5 and ao = 5 [compare Eq. ((3|] at A = 
1.5: In this case, both local spin and local energy densities are 
non-zero during the time-evolution. The inset was reproduced 
from Ref. S 



Therefore, while the analysis of Ref. |28! unveils interest- 
ing details of the time evolution of densities during a Bq 
quench, we maintain that the variance is a useful quan- 
tity to identify candidate parameter sets for spin diffu- 
sion in, e.g., the non-equilibrium regime. Final proof of 
diffusive behavior then needs to be established by either 
demonstrating the validity of the diffusion equation or 
by computing correlation functions, see, e.g., Ref. I55ll6ll . 
For instance, in Ref. [2^ Jesenko and Znidaric analyze 
the steady-state currents in the A > 1 regime at finite 
temperature and obtain diffusive behavior. 



VI. SUMMARY 

We studied the real-time energy dynamics in XXZ 
spin-1/2 chains at zero temperature in two different sce- 
narios. First, we investigated the energy dynamics in the 
absence of spin currents induced by a local sign change 
in the exchange interactions. The spatial variance be- 
haves as 6a'^{t) oc P for all A, consistent with ballis- 
tic dynamics. In the gapless regime, the velocity of the 
fastest excitation present in the dynamics is the group 
velocity Vg of spinous, yet our particular quench also in- 
volves excitations with much smaller velocities resulting 
in expansion velocities Ve < Vg. Furthermore, the bal- 
listic dynamics can be related to properties of energy 
currents. While the total current vanishes in our set-up, 
i.e., (Je) '■= J2i{jf) = 0' tlie current in each half of the 
chain ( Jf'/j) > takes a constant value, after some tran- 
sient dynamics. Therefore, in each half of the system, 
we prepared a state with a conserved global current, al- 
lowing us to make a direct connection to the predictions 
of linear-response theory where the existence of ballistic 
dynamics is directly linked to conservation laws that pro- 
hibit currents from decaying4^ Moreover, we identified 
an observable J^ built from local currents whose expec- 
tation value (o^|;(t)) is time independent if Sa% oc t^ and 
vice versa. This carries over to other types of transport 



as well and, in fact, the analysis of the time-dependence 
of (Je) can be used as an independent means to identify 
ballistic regimes, or to unveil the absence thereof. 

In the second part, we studied the energy dynamics 
induced by quenching a Gaussian magnetic field, with 
two main results. These quenches allow us to access 
the regime of weakly perturbed initial states and in that 
limit, we recover the predictions from Luttingcr Liquid 
theory for the wave-packet dynamics: Their variance sim- 
ply grows as Sa'^ = v^t'^. In the massive phase a very 
interesting phenomenon occurs, since the energy dynam- 
ics is ballistic on time scales over which the spin dy- 
namics behaves diffusively although both are driven by 
the same perturbation. This resembles the picture es- 
tablished from linear-response theory;^ there applied to 
the finite-temperature case, in the non-equilibrium setup 
studied here. While our numerical results cover spin 
chains on real space lattices and initial states far from 
equilibrium, the extension of our work to a finite tem- 
perature of the background will be crucial to tackle the 
most important open questions. 
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Appendix A: Entanglement growth 

Here we want to study the growth of entanglement 
across a junction separating regions in a spin chain with 
ferromagnetic correlations from ones with antiferromag- 
netic ones. 

To that end, we take initial states inspired by Ref. [l9l 
where one half of the system has a positive and the other 
one a negative J. We obtain this configuration as a vari- 
ation of Ji quench choosing: 



J.= 



J i<L/2 

for i = L/2 

-J i > L/2 



(Al) 



in Eq. ^. We then perform the time evolution under the 
antiferromagnetic Hamiltonian [Eq. ([1])]. As a measure of 
the entanglement we calculate the von Neumann entropy 



5'vN = -Tr{pAlnpA) 



(A2) 
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of the reduced density matrix pA = TrBP, where p = 
\ip{t)){ilj{t) and \ip{t)) is the time-evolved wave function, 
for a bipartition in which we cut the chain into two halves 
of length L/2 across the central link. Our results are 
plotted in Fig. [141 We observe that the von-Neumann en- 
tropy grows at most logarithmically (purple dashed line) , 
in agreement with Ref. [2l|. The overall largest values of 
SvNit) are found at the critical point A = 1 (red squares). 
This behavior is very similar to the observations made in 
Ref. Il9| for spin dynamics starting from a state with all 
spins pointing up(down) in the left(right) half. 




10 15 20 

time tj 



FIG. 14: (color online) Time dependence of the von-Neumann 
Entropy Svn for a bipartition that cuts the system across 
the central bond during the time evolution starting from a 
ferromagnetic region coupled to an antiferromagnetic one at 
A = 0.5, 1,1-5 
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